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ABSTRACT 



Context. Although the theoretical understanding of the nonlinear gravitational clustering has greatly advanced in the last decades, in particular 
by the outstanding improvement on numerical N-body simulations, the physics behind this process is not fully elucidated. 
Aims. The main goal of this work is the study of the possibility of a turbulent-like physical process in the formation of structures, galaxies and 
clusters of galaxies, by the action of gravity alone. 

Methods. We use simulation data from the Virgo Consortium, in ten redshift snapshots (from to 10). From this we identify galaxy-sized and 
cluster-sized dark matter haloes, by using a FoF algorithm and applying a boundedness criteria, and study the gravitational potential energy 
spectra. 

Results. We find that the galaxy-sized haloes energy spectrum follows closely a Kolmogorov's power law, similar to the behaviour of dynami- 
cally turbulent processes in fluids. 

Conclusions. This means that the gravitational clustering of dark matter may admit a turbulent-like representation. 

Key words, gravitational clustering - turbulence - numerical N-body simulations - galaxies: structure and evolution - methods: statistical 



1. Introduction 

In a general framework, the standard model of structure for- 
mati on is based o n the gravitational Jeans instability criterion 
(e.g. lLinded[l997l) which neglects viscous, turbulence and dif- 
fusion effects. In this model, the gravity acts to amplify the tiny 
random fluctuations in the primordial density field to the sig- 
nificantly overdense structures we see today (galaxies, groups, 
clusters and superclusters). As long as the density contrast of 
the fluctuations remains small (S <K 1) - in early times for small 
spatial scales (like the ones of galaxies) and until the present 
time for the largest scale structures - the amplitude of these 
fluctuations grow s linearly. In this regi me all scales evolve in- 
dependently (e.g. lSahni & Colesl 1995 ). and the random-phase 
nature of the initial density field is preserved. Analytically, the 
equations of motion of the particles that compose the "cosmic 
fluid" can be successfully solved by means of the perturba- 
tion theory. It turns out that, in the hierarchical scenario, fluc- 
tuations of increasing spatial scale successively enter a non- 
linear regime, collapse and evolve to the virial equilibrium. 
Conversely to the case of the linear regime, there is no gen- 



eral exact ana lytic solution for the evolution of structures on 
this stage (e.g. Padmanabhan 2006). 



* Part of the calculations made use of the NEC SX4/8A supercom- 
puter at the CeNaPAD Ambiental, CPTEC/INPE, Brasil 



The search for a theoretical understanding of the nonlin- 
ear gravitational clustering has greatly advanced in two distinct 
directions: by constructing analytical approximation schemes 
and by evolving numerical N-body simulations. While the at- 
tempts on the first approach are limited to the domains of valid- 
ity of the approximations, the ones on the last reproduce quite 
well most of the observed properties of the actual structures 
(see, for instance. ISpringel et al. Il2005l) . but do not clarify com- 
pletely our understanding of the physics behind it. 

On the analytical front, the most appealing proposals have 
been the ones based on a hydrodynamical approach. The origi- 
nal and simplest mod el of this type is th e one called Zel'dovich 
Approximation (ZA, ZeFdovichl 1970h . which considers that 
the velocity of the fluid elements can be expressed at any time 
in terms of the initial gravitational potential, that is, they move 
in the field without modifying it. Remarkably this model pre- 
dicts the walls (pancakes), filaments, clumps and voids that 
characterize the large scale structure. Nevertheless, the lack 
of self-gravity leads to the dissolution of the structures af- 
ter the first crossing of the particles. Further progress was 
achieved with a series of approximations referred to as ad- 
hesion models (Gurbatov et al. 1989; Shandarin & Zel'dovichl 
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19891: iRibeiro & de Fariatf 2005). Mathematically, these models 



are based on the three-dimensional Burgers' equation dBurgers 



1940) which introduces an ad hoc artificial viscosity that slows 
down the particles as they approach the pancakes and fila- 
ments. The alternative schemes are mostly extensions of the 
ZA, among which are the more recently proposed Lagrangian 
approximations (see, e.g., JSahni&Coleslll995l: iMakler et al.1 
200ltlBernardeau et alJl2002 . and references therein). 

One possibility that was not much explored is that the for- 
mation of structures is more closely related to hydrodynamics, 
for example, by the existence of a "turbulent-like" behaviour 
of the dark matter clustering process. The idea comes from the 
characteristic pattern of the large scale structure that resembles 
the turbulent pattern of fluids, although collisionless dark mat- 
ter is the responsible for such pattern (baryonic matter effects 
take place mostly on smaller scales, while dark energy is ex- 
pected to be acting more recently). 

In this paper we use a N-body simulation from the Virgo 
Consortium to to look for signatures of a turbulent-like process 
in the formation of structures. For this sake we detect systems 
of different scales, from galactic size subhaloes to supercluster 
size super-haloes, study some of their evolutionary properties, 
and evaluate the gravitational energy of bound systems. 



2. The Virgo N-body Simulation 

Our data came from the intermediate scale cosmological N- 
body sim ulations run by the V irgo Consortium^ which are re- 
ported in [Jenkins et al] (11998b . Below we briefly describe the 
main steps and characteristics o f these simulations (for an ex- 
tensive review see, for instance, iBertschingerll 19981 and refer- 
ences therein). 

As usual, the fundamental properties of the N-body simula- 
tions of structure formation in the Universe are defined by the 
background cosmological model and the initial perturbations 
imposed on this background. Four versions of the cold dark 
matter model are available from the Virgo Consortium project, 
from which we chose the ACDM one, with cosmological pa- 
rameters {Q. M , O a , h, o-gfl = (0.3, 0.7, 0.7, 0.9). Among the ad- 
ditional (numerical) parameters that characterize the simulation 
are the simulation box of side 239.5 lr l Mpc and the number of 
particles, 256 3 , with individual masses of 6.86 x 10 10 /; -1 Mq. 
Comoving spatial coordi nates and periodic boundary condi- 



tions were also used (see Jenkin s et al]|1998l for details). Ten 



snapshots are available, atredshifts: 10.0, 5.0, 3.0, 2.0, 1.5, 1.0, 
0.5,0.3,0.1 and 0.0. 

The initial conditions are set by assigning a position and 
velocity to each particle using an algorithm imposing pertur- 
bations on an initially uniform st ate represente d by a "glass" 
distribution of particl es (see, e.g., White 1996). A convenient 
and efficient method dEfstathiou et al. 1985b for converting the 



early density fluctuations into th ese perturbations can be de- 
rived from the Zel'd ovich (1970) solution: 



x = r/a = q + D p 



dr 

dt 



Hr 



dD 



(D 



where x is the comoving Eulerian coordinate of a particle, r 
is its proper Eulerian coordinate, a(t) is the expansion factor 
(which satisfies the Friedmann equation), q is the Lagrangian 
coordinate of the particle denoting its initial unperturbed po- 
sition, D(t) is the growth factor of the linear growing mode, 
p(q) is the displacement field, v is the peculiar velocity of the 
particle and, finally, H{t) = dlna/dt, is the Hubble parame- 
ter. The first term of the displacement equation describes the 
cosmological expansion, while the second denotes the pertur- 
bation. The initial velocity field is defined by the initial density 
fluctuations: 



V.p = -d/D 



(2) 



where S(x, f) = [p - p]/p, p(x, t) is the local mass density and 
p(t) is the spatial mean density. 

Since cosmological N-body simulations are intended to fol- 
low the evolution of the dark matter from the linear regime 
into the deeply nonlinear regime, they must begin with the ex- 
pected distribution of matter, at some early epoch, produced by 
the linear evolution of primordial fluctuations. If these fluctu- 
ations are Gaussian, they can be fully described by the den- 
sity fluctuation power spectrum (PS), P(k, t) = (|(5j:(OI 2 >- This 
PS can be written as the product of an initial power law, pro- 
duced by the process which generated the fluctuations (typ- 
ically the amplification of quantum fluctuations during infla- 
tion), a transfer function T(k) representing the linear evolution 
of each mode through the early expansion history (matter dom- 
ination, baryon-photon decoupling, etc.), and the growth factor 
D(f) governing the later linear evolution after decoupling: 



P(k,t)=Ak"T(k) 2 D(t) 2 



(3) 



where A is a normalization factor encoding the overall ampli- 
tude of the initial fluctuations (generally determined from erg). 

For the Virgo simulations, a Harrison-Zel'dovich n — 1 
spectral index was used, and also the approximation given by 
Bond & Efstathioul 0984) for the transfer function: 



7Xk) = {l + [a4 + (b4) 3/2 + (c<7) z ] 



-,2-ivi-l/v 



(4) 



1 available from the web-page 

http : / /www . mpa-garching . mpg . de/Virgo/data_download . html 

2 a M = P/Pc = ZnGplQH 2 ); 
n A = Ac 2 /(3// 2 ); 

h = HalWOkms- 1 Mpc' 1 ; and 

cr 8 is the rms mass fluctuation within a top-hat radius of 8h~ l Mpc. 



where q = k/T, a = 6.4/T 1 Mpc, b = 3 ft -1 Mpc, c = 1.7ft" 1 
Mpc and v = 1.13. So, the shape of the PS is determined by 
the parameter F = Qrj/j. The amplitude A, on the other hand, 
was obtained from the observed abundance of galaxy clusters 
dJenkins et alJll998h . 

Given the initial conditions, that is, the initial position and 
velocity vector of each particle, the simulation is run in about 
one thousand timesteps from z = 30 to z = 0. Newton's law 
and Poisson's equation, in commoving coordinates, are used to 
evolve forward in time the particle trajectories: 



-V<p(x), 



<P(x) 



N 



j-f [(x-xtf + e 1 ] 



(5) 
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The parameter e is a softening factor introduced in order 
to prevent the formation of unphysical binaries and large an- 
gle particle scatterings, and to ensure that the two-body re- 
laxation time i s sufficiently large to guarant ee the collision- 
less behaviour (Efstath iou & Eastwood! Il98ll) . The time inte- 
gration is performed u sing a second order leapfrog scheme 
(Efstathiou et al .] Il985h . At each timestep the gravitational 
force of every particle, generated by all others, must be calcu- 
lated. This is a critically time consuming task if done by direct 
summation, so that simplifying parallel algorithms may be used 
to efficiently compute it. The Virgo Consortium used a parallel 
adaptive P 3 M algorithm (jCouchman. Thomas & Pearce 1995 : 



Pear ce & Couchmanll 19971) . which computes long range grav- 



itational forces by smoothing the mass distribution in a mesh 
(particle-mesh, PM), short range gravitational forces by direct 
summation (particle-particle, PP), and applying high resolution 
refinements around strong clustered regions. 

3. Identification of dark matter structures 

There are many methods for detecting dark matter haloes 
in N-body simulations. The most simple, used and tested 
ones are the p ercolation algorit h m (or "Friends-of-Fr i ends" , 
dubbed FoF iHuchra & Gellerl Il982t lEinasto et al.1 1 1984 ; 
)avisetal] 1 19851) the "Spherical Overdensity" method 
Lacev & Cold 1 1994 and the "DENMAX" sliding scheme 



Gelb & Bertschinger1l994 



The FoF algorithm groups together particles that have pair 
separations smaller than a chosen linking length, I. This linking 
length is frequently referred to as b times the mean interparticle 
separation. The resulting "groups" are bounded by a surface of 
approximately constant density: 



1 



(4/3);rf 3 n 2nfi 



2nb 3 



I 

2bi 



Assuming that the density profile of these groups can be ap- 
proximated by isothermal spheres, the average density contrast 
internal to this surface is given by about 3 times the surface 
density. The larger is I (or b), the lower is the density con- 
trast and the higher is the number of particles linked to the 
groups. In general the value of b is chosen to give a mean 
overdensity close to the value expected for a virialized ob- 
ject in the framework of the spherical collapse model, ~179. 
This gives b = 0.2, which see ms to be an appropriate value for 



detecting cluster-sized haloes (ILacey & Cold 1 1994: Eke et al 



1996; Governato et al . 1999). The main advantages of the FoF 



method are the simplicity (only one free parameter), the re- 
produptibility (a unique group catalog for any chosen value of 
€) and the capability of detecting haloes of any shape. 

The other approach, SO algorithm, identifies density peaks 
and puts spheres around them, increasing the radius of the 
spheres until the average density contrast reaches a chosen 
value. It was proposed to overcome the problem of FoF that 
may accidently link two distinct lumps if they are connected 
by a low density bridge of particles. However, it tends to loose 
the outer portions of el lipsoidal haloes due to the assumption 
of spherical symmetry dLacev & Coldll994 . 



The DENMAX algorithm, on the other hand, "moves" par- 
ticles along local density gradients toward density maxima, 
separating halo candidates by three-dimensional density val- 
leys. The improvement of this method is the application of a 
self-boundedness check by eliminating particles with positive 
total energy. Unfortunately, this scheme has the drawback that 
the results depend on the lev el of smoothing of the density field 
dGelb & Bertschingerlll994 . 

Nevertheless, the three methods abo v e are found to 
give similar resu lts dLacev & Cold Il994t lEke et al. 1996; 



lAudit et all fl998h . The main limitation of these algorithms 
for detecting structures on N-body simulations appears for 
small masses (like the ones of galaxies): they do not re- 
solve properly "sub-haloes" embedded in host (cluster-sized) 
haloes with the usual parameters. Improved versions of 
FoF, S O and DENMAX - respectively the hierarchical FoF 
rHFoF.lKlvpin et alJl999h . the B ound Density Maxima ( BDM, 
Klypin et al J 19991) and the HOP dEisenstein & Hut! 19981) - and 
also t wo step algorithms - like the SU BFIND dSpringel et alJ 



200 lb and the Physically Self-Bound dKim & Parkl 120061) - 



have been proposed to better identify sub-haloes in crowded 
regions. 

Here we use a different approach. It works, in some sense, 
like the HFoF algorithm. Since the FoF with b = 0.2 was shown 
to successfully detect cluster-sized haloes, we use this scheme. 
A cut in a minimum number of particles is used to eliminate 
haloes with masses smaller than the ones expected for cluster- 
sized haloes (see table 1). For galaxy-sized haloes (including 
sub-haloes) we use FoF with a smaller linking parameter (b = 
0.1). This results in a catalog of haloes smaller and denser than 
the cluster-sized haloes. Also, the equivalent mean density con- 
trast (~1 300) is close to the one expected for the Milky-Way 
if one considers the mass estimate of 1.5 x 10 12 hr x Mp) and 



the radius of 150 h 1 kpc (e .g. ISmith et alJl2007t iDehnen et al 



(6) 120061: ISakamoto et alJl2003b . So, this is appropriate for identi 



fying sub-haloes and also for isolated galaxy-sized haloes. cD 
galaxies are the most massive galaxies and, since they are tran- 
sition objects (their halo may be indistinguishable from their 
parent cluster halo), they are detected as overmassive galaxy- 
sized haloes. They also coincide with their parent cluster po- 
sition. For some analyses we excluded these special galaxies 
by putting a maximum mass limit. We also tentatively identify 
supercluster-sized haloes by relaxing the linking parameter to 
b — 1.0 (density contrast about 1.5). Superclusters are not viri- 
alized systems, but structures that probably have just detached 
from the Hubble flow. The volume of the simulation is small 
for accurate statistics of superclusters, but we can select some 
typical cases for further analysis (to be published elsewhere). 

We also applied a boundedness check for galaxy-sized 
haloes, evaluated in a statistical way. FigureQ]shows the distri- 
bution of galaxy-sized halo masses (by the number of grouped 
particles) with their respective velocity dispersions. From the 
Virial Theorem we expect that, for relaxed bound structures, 
these two parameters are correlated, that is, the velocity disper- 
sion is proportional to the square root of the mass. The ridge 
of points in the lower part of the figure follows closely this 
relation, as can be seen by the fit (solid line). So, we consid- 
ered as bound the haloes below the dashed line (about 2.5 rms 
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Table 1. Halo identification and selection parameters (/ is 
the mean density contrast). 



Halo size 


b 


I 


/ 


Mass 






(Mpc) 




(N m ) 


(M Q ) 


Galaxies 


0.1 


0.10 


1300 


2-150 


> 1 x 10" 


Groups/Clusters 


0.2 


0.18 


180 


> 100 


> 7 x 10 12 


Superclusters 


1.0 


0.92 


1.5 


> 15 000 


> 1 x 10 15 



:;::(; 



800 



400 



300 




150 



part 



Fig. 1. Distribution of galaxy-sized halo masses (N part ) with 
their respective velocity dispersion (<x/ M / ), for redshift 0. The 
solid line is the fit to the points in the main ridge (a power law 
Chaio K N parl ^ 2 is used), while the dashed one represents 2.5 
rms of dispersion above this fit. 

above the fit). For z = 0.0, 50.7% of the particles remained 
isolated after the application of the FoF algorithm to search 
for galaxy-sized haloes, and 10.6% were excluded as unbound 
haloes. These isolated particles and the particles of excluded 
haloes are probably associated to galaxies smaller then our de- 
tection limit. A similar behaviour was found for the other red- 
shift snapshots. 

4. Correlation Functions 

In this section we discuss the clustering properties of the 
galaxy-sized and cluster-sized haloes, found by our percola- 
tion scheme, in order to check the consistency of our object 
catalogs. We measure the t wo-po int (auto)correlation function 
using the lLandy & Szalayl (119931) estimator: 



= -EBWD(n R ln D f - 2DR(n R /n D ) + RR] 



(7) 




r (h 1 Mpc) 
Fig. 2. Correlation functions for dark matter particles. 



and n D are the mean number densities in data and random 
catalogs. We used n R — n n in all cases. We also applied the 
Davis & Peebles! d 1983b and lHamiltonl d 1993b estimators to our 
data and found very similar results on the sampled ranges. 

Figure [2] shows the two point correlation functions for the 
dark matter particles in the simulation (DMCF), for redshifts 
between 0.0 and 5.0. The results for redshi fts 0.0, 0.3, 0.5 and 
1 .0 are very close to the ones presented by lJenkins et al.1 (11998) 
on their third top panel of figure 5, as expected. The z — 0.1 
correlation function has a slightly different behaviour on scales 
larger than 1 hr l Mpc. It is a well established fact that the dark 
matter (mass) correlation function has a very complex shape 
and canno t be represented by a single power law in a significant 
range (e.g. Jenkins et aljl998t Pearce et al.ll999l Benson et al 



where DD, DR and RR are the pair counts, respectively in the 
data-data, data-random and random-random catalogs, and n R 



2000l: ISpringel et al.ll2005l) . Another known property that can 
be clearly seen is that the amplitude of the DMCF evolves con- 
siderably with time. Also its slope changes systematically with 
z. 

The correlation functions for galaxy-sized haloes (GHCF) 
are presented on Fig. [3] for the same range of redshifts. Unlike 
the DMCF, it shows a shape that resembles closely a power law 
in all redshifts and over a large range of pair separations (here 
we display the range adequately sampled by our data, 0.5 < 
r < 20 Mpc). This behaviour has been pointed as r eally no 
more than a coincidence (see, e.g. Sprin gel et al ]|2005b . but, as 
we shall see in Sect. 5, this might be related to the emergence of 
a turbulent-like process in the halo clustering. On small scales 
(less than a few Mpc) and low redshifts, the GHCF has less 
power than the DMCF, and the GHCF is said to be "antibiased". 
On larger scales, on the other hand, the GHCF function show 
similar amplitudes to the DMCF for lower redshifts, and so it 
is "unbiased". Furthermore, the GHCF evolves very little with 
redshift (its amplitude grows slightly with time). 
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r (h Mpc) 



Fig. 3. Correlation functions for galaxy-sized haloes, 
types are the same as on Fig. 2. 



Line 




r (h 1 Mpc) 

Fig. 4. Correlation functions for cluster-sized haloes. Line 
types are the same as on Fig. 2. 



Figure|4]displays the correlation functions for cluster-sized 
haloes (CHCF). One can readily see that, up to redshifts z ~ 1, 
it evolves even slower than the GHCF At higher redshifts, the 
CHCF begins to grow. This is a consequence of the correlation 
between the correlation function of clusters and their masses 
or richness: at higher redshifts, only the m ost massive clusters 
are detected and they are more cl ustere d dMo & Whitd fl996; 
Bahcall et al]|2004UYounger et alj|2005h . 
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Fig. 5. Evolution of the parameters amplitude [panel (a)] and 
slope [panel (b)] for the power law fit to the correlation func- 
tions of galaxy-sized haloes and cluster-sized haloes. The dot- 
ted lines are the linear fits to the data. 

In order to better characterize the evolution of the correla- 
tion strength (ro) and slope (y) with redshift we fitted power 
laws to the correlation functions of galaxy-sized and cluster- 
sized haloes 



(8) 



respectively in the ranges 0.6-18 h' 1 Mpc and 2-18 h Mpc. 
The results can be seen on Fig. [5] Panel (a) shows the evo- 
lution of ro, which grows slightly with time for both galaxy- 
sized haloes and cluster-sized haloes (more for the former), ex- 
cept for the high redshift points for which the biased sampling 
towards more massive objects, as cited above, begins to dom- 
inate. Panel (b) displays the y evolution. There is almost no 
change in the slope, except for probable random noise that also 
grows with redshift due to the decrease in the number of ob- 
jects. The mean slope for galaxy-sized haloes and cluster-sized 
haloes are, respectively, 1.53 and 1.59, with standard devia- 
tions of 0.03 and 0.09. This values are smaller than the classical 
1.8 slope, but consistent with recent observational results, espe 



dally for galaxies (e.g. Padilla & Baugh 2003j lHawkins et al 
2003h . 



5. Turbulence-like Energy Spectra 

Now we test the hypothesis of the existence of turbulent-like 
processes in the formation of structures, searching for their 
possible signatures. The main characteristic of 3D turbulence 
is the existence of a kinetic energy cascade. According to this 
phenomenology, kinetic energy is injected in the fluid at large 
scales by an external mechanical forcing. Large scale eddies 



6 



C. A. Caretta et al.: Turbulence on Formation of Galaxies 







1 


im 


1 1 1 1 III 

-^5/3 

■ ■ . . * • . - 




1 1 

Clusters 




100 










* » . ■ . • » " • 






2 










°°°°°o„_ 


• 

• 


m 

• ■ 

• 


o 






□ □ 


□ □ 


° 




















r-T n q ^ ° ° * O 








10 










p" ; ; ; ; ■ ■ 

i i i 

i 


o 

i 





0.1 

k (Mpc -1 ) 



0.1 




k (Mpc -1 ) 



Fig. 6. Gravitational energy spectra for cluster-sized haloes (above panel) and galaxy-sized haloes (below panel), for the sampled 
redshifts. A reference -5/3 power law is plotted as a solid line. 



are then deformed and stretched by the fluid dynamics, break- 
ing into smaller eddies. This process is repeated, through a 
hierarchy of smaller and smaller eddies, until the kinetic en- 
ergy is finally dissipated into heat by the viscosity of the fluid 
dFrischl 1 1 995D . Within a certain range of intermediate scales, 
the distribution of energy among turbulence vortices behaves 
like a power law with exponent -5/3. This is the celebrated 
Kolmogorov's energy spectrum, and represents the key signa- 
ture of the "top-down" turbulent cascade scenario. Naturally, 
in an unbounded medium filled with a collisionless set of par- 
ticles, the emergence of a turbulent-like behaviour would have 
a different physical origin than in the standard hydrodynamic 
turbulence. Nevertheless, the effects of particle trajectory cross- 



ings and virialization, as well as the damping effect of the uni- 
versal expansion, could lead to similar phenomena as observed 
in standard turbulence. For example, the cosmological dynam- 
ics of structure formation could exhibit the power law scaling. 
This will be the focus of this section. In other words, we will 
investigate the energy distribution of the cosmic structures dis- 
cussed above (i.e. galaxies and clusters) to look for the signal 
of such a turbulent-like behaviour in the form of a power law 
scaling. 

Since we are not dealing with an usual fluid (dark matter is 
acollisional and, so, not viscous) and we have no energy "injec- 
tion" in the strict sense, our assumption must be that the energy 
that produces the turbulent-like behaviour is the gravitational 
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potential energy. Thus, with the lists of identified galaxy-sized 
and cluster-sized haloes we obtained the energy spectra of the 
gravitational potential. This was done by estimating, for each 
halo and in concentric shells, the gravitational energy due to all 
other haloes: 



the use of appropriate measurement techniques (iRamos et al 
2002r.lWuensche et alJEooitlAndrade et alJl2006h . 



1 v-i rrij 

u J=? Gm jL-r. 



(9) 



where ry is the distance between the haloes i and /'. Then we 
estimate the mean energy (for all haloes) of each shell and 
build the energy spectra with its cumulative value. The re- 
sults of this calculation are presented on Fig. [6] panel (a) for 
cluster-sized haloes and panel (b) for galaxy-sized haloes. The 
stochastic error estimations (1/ y/N) are plotted only for the 
first and last redshifts to make the viewing of the figure more 
clear. The -5/3 power law is represented by a solid line. One 
can easily notice that the energy spectra for the cluster-sized 
haloes, although following a power law, does not exhibit the 
Kolmogorov's index. On the other hand, the one for galaxy- 
sized haloes follows closely the Ar 5 ^ 3 scaling in the wave num- 
ber range 0.02 to 0.07 (15 to 50 hr x Mpc). 

It should be noted that this result differs from the one ob- 
tained for highly compressible turbulence in molecular clouds 



( Krit suk et al.l 2007; Padoa n et alJ 120071) . Three-dimensional 



simulations of supersonic Euler turbulence show a velocity 
power spectrum that gets steeper as the Mach number in- 
creases, reaching the Burgers slope of -2 asymptotically. 
Kolmogorov's -5/3 scaling is only recovered by mixing the 
velocity and density statistics, through the computation of 
density-weighted velocity spectra. The very different context 
of the present paper (namely, gravitational clustering of large- 
scale structures, N-body simulation, potential energy spectra) 
may explain this discrepancy. 

To interpret this behaviour we point out that the turbulent- 
like spectra shown in Fig. [6] have only the inertial regime, 
without both the typical spectral energy-input and dissipa- 
tion signatures, which are usually fou nd in standard turbulent 
system s. As originally suggested by Shandarin & Zel'dovichl 
( 1989b . a simple turbulent scenario based on the ZA is com- 
patible with a top-down structure formation, where a large 
structure (pancake), corresponding to the integral scale, forms 
first then fragments into smaller objects respecting a cascade 
process. However, the observations and simulations have, for 
many years, favored a scenario closer to the hierarchical one. 
In such a way, since the gravitational potential energy is purely 
attractive, its variation as a function of the "wave number", de- 
fined from the distances among haloes, could be interpreted 
as a bottom-up cascade process. Physically, because of this at- 
tractive nature of gravity, there is always an intermittent in- 
trinsic instability that amplifies any small homogeneity devi- 
ation, as does the equivalent hydrodynamic turbulence. Thus, 
interpreting the results of Fig. [6] as the outcome of a kind 
of "gravitational turbulence" means that the expansion of the 
universe in the largest scales is driving an inverse fragmen- 
tation process to compensate its deviation from homogeneity 
in local scales along its evolution. Furthermore, a turbulent- 
like behaviour may be a typical and essential multi-scaling 
cosmological feature, requiring, for its robust characterization, 



6. Concluding Remarks 

We identified galaxy-sized and cluster-sized haloes in an in- 
termediate scale ACDM cosmological N-body simulation by 
the Virgo Consortium. For these objects we calculated a grav- 
itational potential energy spectra and found that the one for 
galaxy-sized haloes may be well described by a -5/3 power 
law in the range from 15 to 50 h Mpc. It is interesting to 
note that the same scaling has been found in the observational 
front, for elliptical g alaxies, between potential energy and mass 
(Marquez et alj|2001 ). We speculate that this scaling of the en- 
ergy spectrum, shown in Fig. [6] could be a signatur e of the 
Kolmogorov's universality assumption ( Frischl 19951) . where 
the system dynamical structural behaviour is uniquely and uni- 
versally determined by the scaling and its associated mean en- 
ergy rate. From our results, the formation of structures in the 
nonlinear regime seems to be driven by a turbulent-like mecha- 
nism, associated with irregular fluctuations in an unstable grav- 
itational field, characterized by combinations of small-scale ed- 
dies and larger flow-like structures. 

It should be noted that a power-law spectrum by itself does 
not mean the existence of turbulence. However, the emergence 
of a (approximately) scale invariant hierarchy of structures as 
the result of the gravitational clustering process, together with 
the results of Fig. [6] are consistent with the turbulence-like sce- 
nario suggested here. Also, the space-time patterns seen in the 
simulated data can be interpreted as intermittent nonlinear fil- 
aments in a turbulent-like "gravitational" fluid (see, for phe- 
nomenol ogical compari son purposes, the turbulent fluid simu- 
lated by lShe et al.ll 199 lr) . Because of its importance in our ap- 



proach, it is worth mentioning that there is presently no closed 
theoretical description of turbulence 



ire is p 
(e.gj 



Sreenivasan 1995 



Velho et alfcOOlHRamos et al.ll200lh . Detailed studies, includ- 
ing analyses of higher resolution cosmological simulations and 
theoretical advances on the models for the nonlinear regime 
and turbulence itself, are needed to advance on interpreting 
the apparent turbulent-like behaviour of the structure formation 
processes presented here. 
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